
##### Replication script for the article: Pathways for overcoming commitment problems in the transition from intra-state conflict to peace: a fuzzy set qualitative comparative analysis of 34 comprehensive peace agreements ###

### Ibrahim Kumek ###




### COMMANDS BEFORE THE ANALYTIC MOMENT ###

# Install packages
install.packages(c("psych", "QCA", "SetMethods"), dependencies = TRUE)


# Load packages
library(psych); library(QCA); library(SetMethods)


# Set working directory
setwd("my_workingdirectory")


# Read a CSV formatted dataset using a comma as decimal point and a semicolon as field separator
mydata_cpa_id <- read.csv("mydata_cpa_id.csv", row.names=1, header = TRUE, sep = ";",  dec = ",")


# Save a CSV formatted dataset using a comma as decimal point and a semicolon as field separator
write.csv2(mydata_cpa_id, "mydata_cpa_id.csv") 


# View the dataset 
View(mydata_cpa_id)


# View the first six rows of the dataset
head(mydata_cpa_id)


# Check how many variables there are in the dataset
length(mydata_cpa_id)


# Get the names of the variables
colnames(mydata_cpa_id)


# Get the names of the cases
rownames(mydata_cpa_id)


# Descriptive statistics for the dataset
describe(mydata_cpa_id)




### COMMANDS FOR THE CALIBRATION PROCESS ###

# Looking up the raw values
mydata_cpa_id$power_sharing

mydata_cpa_id$transitional_justice

mydata_cpa_id$transitional_reform

mydata_cpa_id$international_assistance

mydata_cpa_id$implementation_process

mydata_cpa_id$conflict_years_15


# Calibrating fuzzy sets with direct method 
mydata_cpa_id$PS <- calibrate(mydata_cpa_id$power_sharing, type = "fuzzy", thresholds = c(0.05, 1.5, 2.95), logistic =TRUE)
mydata_cpa_id$PS

mydata_cpa_id$TJ <- calibrate(mydata_cpa_id$transitional_justice, type = "fuzzy", thresholds = c(0.05, 9.5, 11.95), logistic =TRUE)
mydata_cpa_id$TJ

mydata_cpa_id$TR <- calibrate(mydata_cpa_id$transitional_reform, type = "fuzzy", thresholds = c(0.05, 12.5, 15.95), logistic =TRUE)
mydata_cpa_id$TR

mydata_cpa_id$IA <- calibrate(mydata_cpa_id$international_assistance, type = "fuzzy", thresholds = c(0.05, 1.5, 2.95), logistic =TRUE)
mydata_cpa_id$IA

mydata_cpa_id$IP <- calibrate(mydata_cpa_id$implementation_process, type = "fuzzy", thresholds = c(0.05, 79.5, 99.95), logistic =TRUE)
mydata_cpa_id$IP

mydata_cpa_id$SCNPT <- calibrate(mydata_cpa_id$conflict_years_15, type = "fuzzy", thresholds = c(14.95, 0.5, 0.05), logistic =TRUE)
mydata_cpa_id$SCNPT


# Comparing raw and calibrated scores 
mydata_cpa_id[c("power_sharing", "PS")]

mydata_cpa_id[c("transitional_justice", "TJ")]

mydata_cpa_id[c("transitional_reform", "TR")]

mydata_cpa_id[c("international_assistance", "IA")]

mydata_cpa_id[c("implementation_process", "IP")]

mydata_cpa_id[c("conflict_years_15", "SCNPT")]


# Visualizing with plot 
plot(mydata_cpa_id$power_sharing, mydata_cpa_id$PS, main = "Calibration patterns of PS", xlab = "Raw score", ylab = "Calibrated score")
abline(h = 0.5, v = 1.5)

plot(mydata_cpa_id$transitional_justice, mydata_cpa_id$TJ, main = "Calibration patterns of TJ", xlab = "Raw score", ylab = "Calibrated score")
abline(h = 0.5, v = 9.5)

plot(mydata_cpa_id$transitional_reform, mydata_cpa_id$TR, main = "Calibration patterns of TR", xlab = "Raw score", ylab = "Calibrated score")
abline(h = 0.5, v = 12.5)

plot(mydata_cpa_id$international_assistance, mydata_cpa_id$IA, main = "Calibration patterns of IA", xlab = "Raw score", ylab = "Calibrated score")
abline(h = 0.5, v = 1.5)

plot(mydata_cpa_id$implementation_process, mydata_cpa_id$IP, main = "Calibration patterns of IP", xlab = "Raw score", ylab = "Calibrated score")
abline(h = 0.5, v = 79.5)

plot(mydata_cpa_id$conflict_years_15, mydata_cpa_id$SCNPT, main = "Calibration patterns of SCNPT", xlab = "Raw score", ylab = "Calibrated score")
abline(h = 0.5, v = 0.5)


# Visualizing with histogram, combining two graphs
par(mfrow=c(1, 2))
hist(mydata_cpa_id$power_sharing)
hist(mydata_cpa_id$PS)
par(mfrow=c(1, 1))

par(mfrow=c(1, 2))
hist(mydata_cpa_id$transitional_justice)
hist(mydata_cpa_id$TJ)
par(mfrow=c(1, 1))

par(mfrow=c(1, 2))
hist(mydata_cpa_id$transitional_reform)
hist(mydata_cpa_id$TR)
par(mfrow=c(1, 1))

par(mfrow=c(1, 2))
hist(mydata_cpa_id$international_assistance)
hist(mydata_cpa_id$IA)
par(mfrow=c(1, 1))

par(mfrow=c(1, 2))
hist(mydata_cpa_id$implementation_process)
hist(mydata_cpa_id$IP)
par(mfrow=c(1, 1))

par(mfrow=c(1, 2))
hist(mydata_cpa_id$conflict_years_15)
hist(mydata_cpa_id$SCNPT)
par(mfrow=c(1, 1))


# Select relevant sets, i.e. the calibrated sets from the raw dataset
cpaid_calib <- mydata_cpa_id[ , c("PS", "TJ", "TR", "IA", "IP", "SCNPT")] 


# Saving calibrated data 
write.csv2(cpaid_calib, "cpaid_calib.csv")


# Skewness checks 
skew.check(cpaid_calib) 


# Case names that are more inside than outside of the set
rownames(subset(cpaid_calib, PS > 0.5)) 

rownames(subset(cpaid_calib, TJ > 0.5)) 

rownames(subset(cpaid_calib, TR > 0.5)) 

rownames(subset(cpaid_calib, IA > 0.5)) 

rownames(subset(cpaid_calib, IP > 0.5)) 

rownames(subset(cpaid_calib, SCNPT > 0.5)) 


# Case names that are more outside than inside of the set
rownames(subset(cpaid_calib, PS < 0.5))

rownames(subset(cpaid_calib, TJ < 0.5)) 

rownames(subset(cpaid_calib, TR < 0.5)) 

rownames(subset(cpaid_calib, IA < 0.5)) 

rownames(subset(cpaid_calib, IP < 0.5)) 

rownames(subset(cpaid_calib, SCNPT < 0.5))


# Checking for cases on the crossover point
ambig.cases(cpaid_calib$PS)

ambig.cases(cpaid_calib$TJ)

ambig.cases(cpaid_calib$TR)

ambig.cases(cpaid_calib$IA)

ambig.cases(cpaid_calib$IP)

ambig.cases(cpaid_calib$SCNPT)




### COMMANDS FOR NECESSITY ANALYSIS ###

# Testing necessity of multiple conditions (also calculates necessity of negated conditions) 
QCAfit(cpaid_calib[, c("PS", "TJ", "TR", "IA", "IP")], cpaid_calib$SCNPT, necessity = TRUE)


# Necessity XY plot of multiple conditions
xy.plot(x = "PS" , y = "SCNPT" , data = cpaid_calib, necessity = TRUE , jitter = TRUE , xlab = "PS" , ylab = "SCNPT" , crisp = FALSE)

xy.plot(x = "TJ" , y = "SCNPT" , data = cpaid_calib, necessity = TRUE , jitter = TRUE , xlab = "TJ" , ylab = "SCNPT" , crisp = FALSE)

xy.plot(x = "TR" , y = "SCNPT" , data = cpaid_calib, necessity = TRUE , jitter = TRUE , xlab = "TR" , ylab = "SCNPT" , crisp = FALSE)

xy.plot(x = "IA" , y = "SCNPT" , data = cpaid_calib, necessity = TRUE , jitter = TRUE , xlab = "IA" , ylab = "SCNPT" , crisp = FALSE)

xy.plot(x = "IP" , y = "SCNPT" , data = cpaid_calib, necessity = TRUE , jitter = TRUE , xlab = "IP" , ylab = "SCNPT" , crisp = FALSE)


# Testing necessity of multiple conditions for the negated outcome (also calculates necessity of negated conditions)
QCAfit(cpaid_calib[, c("PS", "TJ", "TR", "IA", "IP")], 1-cpaid_calib$SCNPT, necessity = TRUE)


# Testing for necessary disjunctions
SUIN_cpa <- superSubset(data = cpaid_calib, outcome = "SCNPT", conditions = c("PS", "TJ", "TR", "IA", "IP"), incl.cut = 0.9, cov.cut=0.6, ron.cut = 0.5)
SUIN_cpa

SUIN_cpa[["incl.cov"]]

SUIN_cpa[["coms"]]


# XY plots for necessary disjunctions
pimplot(data = cpaid_calib, results = SUIN_cpa, outcome = "SCNPT", necessity = TRUE, jitter = TRUE, crisp = FALSE, all_labels = TRUE)




### COMMANDS FOR SUFFICIENCY ANALYSIS ###  

# Constructing a truth table for the positive outcome
cpaid_TT <- truthTable(data = mydata_cpa_id, outcome = "SCNPT", 	
                       conditions = c("PS", "TJ", "TR", "IA", "IP"), 
                       incl.cut = 0.8, pri.cut = 0.7, n.cut = 1, 
                       sort.by = c("OUT", "incl"), complete = TRUE, 
                       show.cases = TRUE)
cpaid_TT

cpaid_TT[["tt"]]  

cpaid_TT[["indexes"]]

cpaid_TT[["cases"]]

cpaid_TT[["DCC"]]

cpaid_TT[["minmat"]]


# Constructing a truth table for the negative outcome
cpaid_TTnot <- truthTable(data = mydata_cpa_id, outcome = "~SCNPT", 	
                          conditions = c("PS", "TJ", "TR", "IA", "IP"), 
                          incl.cut = 0.8, pri.cut = 0.7, n.cut = 1, 
                          sort.by = c("OUT", "incl"), complete = TRUE, 
                          show.cases = TRUE)
cpaid_TTnot

cpaid_TTnot[["tt"]]


# Logical minimization for the positive outcome
sol_out <- minimize(cpaid_TT,
                    details = TRUE)
sol_out

sol_out[["tt"]]

sol_out[["negatives"]]

sol_out[["initials"]]

sol_out[["PIchart"]]

sol_out[["primes"]]

sol_out[["solution"]]

sol_out[["essential"]]

sol_out[["inputcases"]]

sol_out[["pims"]]

sol_out[["IC"]]

sol_out[["numbers"]]


# Logical minimization for the negative outcome
sol_nout <- minimize(input = cpaid_TTnot, details = TRUE)
sol_nout


# XY plot sufficiency for the positive outcome
pimplot(data = mydata_cpa_id, outcome = "SCNPT", results = sol_out, all_labels = TRUE, jitter = TRUE, markers = TRUE)


# Radar charts
QCAradar(results = sol_out, outcome = "SCNPT")


# Conservative solution
sol_outc <- minimize(cpaid_TT, details = TRUE)
sol_outc


# Most parsimonious solution
sol_outp <- minimize(cpaid_TT, details = TRUE, include = "?")
sol_outp

sol_outp[["PIchart"]]

sol_outp[["primes"]]

sol_outp[["solution"]]

sol_outp[["essential"]]

sol_outp[["inputcases"]]

sol_outp[["pims"]]

sol_outp[["IC"]]

sol_outp[["SA"]]


# XY plots for most parsimonious solution formula and sufficient terms
pimplot(data = mydata_cpa_id, outcome = "SCNPT", results = sol_outp, all_labels = TRUE, jitter = TRUE, markers = TRUE)


# Radar charts for most parsimonious solution formula and sufficient terms
QCAradar(results = sol_outp, outcome = "SCNPT")


# Venn diagram for the truth table
library(venn)

venn(cpaid_TT, counts = TRUE)


# Intermediate solution
sol_outi <- minimize(cpaid_TT, details = TRUE, include = "?", dir.exp = "TR, IP", row.dom = TRUE)
sol_outi

sol_outi[["i.sol"]]


# XY plots for intermediate solution formula and sufficient terms
pimplot(data = mydata_cpa_id, outcome = "SCNPT", results = sol_outi, all_labels = TRUE, jitter = TRUE, markers = TRUE)


# Radar charts for intermediate solution formula and sufficient terms
QCAradar(results = sol_outi, outcome = "SCNPT")




### COMMANDS FOR ROBUSTNESS TESTS AND CLUSTER DIAGNOSTICS ###  

# Calculating sensitivity ranges for the calibration of PS condition
rob.calibrange(raw.data = mydata_cpa_id,
               calib.data = cpaid_calib,
               test.cond.raw = "power_sharing",
               test.thresholds = c(0.05, 1.5, 2.95),
               step = 0.5,
               max.runs = 20,
               outcome = "SCNPT",
               conditions = c("PS", "TJ", "TR", "IA", "IP"),
               incl.cut = 0.80,
               n.cut = 1,
               include = "?")


# Calculating sensitivity ranges for the calibration of TJ condition
rob.calibrange(raw.data = mydata_cpa_id,
               calib.data = cpaid_calib,
               test.cond.raw = "transitional_justice",
               test.thresholds = c(0.05, 9.5, 11.95),
               step = 0.5,
               max.runs = 20,
               outcome = "SCNPT",
               conditions = c("PS", "TJ", "TR", "IA", "IP"),
               incl.cut = 0.80,
               n.cut = 1,
               include = "?")


# Calculating sensitivity ranges for the calibration of TR condition
rob.calibrange(raw.data = mydata_cpa_id,
               calib.data = cpaid_calib,
               test.cond.raw = "transitional_reform",
               test.thresholds = c(0.05, 12.5, 15.95),
               step = 0.5,
               max.runs = 20,
               outcome = "SCNPT",
               conditions = c("PS", "TJ", "TR", "IA", "IP"),
               incl.cut = 0.80,
               n.cut = 1,
               include = "?")


# Calculating sensitivity ranges for the calibration of IA condition
rob.calibrange(raw.data = mydata_cpa_id,
               calib.data = cpaid_calib,
               test.cond.raw = "international_assistance",
               test.thresholds = c(0.05, 1.5, 2.95),
               step = 0.5,
               max.runs = 20,
               outcome = "SCNPT",
               conditions = c("PS", "TJ", "TR", "IA", "IP"),
               incl.cut = 0.80,
               n.cut = 1,
               include = "?")


# Calculating sensitivity ranges for the calibration of IP condition
rob.calibrange(raw.data = mydata_cpa_id,
               calib.data = cpaid_calib,
               test.cond.raw = "implementation_process",
               test.thresholds = c(0.05, 79.5, 99.95),
               step = 5,
               max.runs = 20,
               outcome = "SCNPT",
               conditions = c("PS", "TJ", "TR", "IA", "IP"),
               incl.cut = 0.80,
               n.cut = 1,
               include = "?") 


# Calculating sensitivity ranges for the raw consistency threshold
rob.inclrange(data = mydata_cpa_id,
              step = 0.01,
              max.runs = 40,
              outcome = "SCNPT",
              conditions = c("PS", "TJ", "TR", "IA", "IP"),
              incl.cut = 0.8,
              n.cut = 1,
              include = "?")


# Calculating sensitivity ranges for the frequency threshold
rob.ncutrange(data = mydata_cpa_id,
              step = 1,
              max.runs = 20,
              outcome = "SCNPT",
              conditions = c("PS", "TJ", "TR", "IA", "IP"),
              incl.cut = 0.8,
              n.cut = 1,
              include = "?")


# Calculate robustness parameters
IS <- minimize(data = mydata_cpa_id,
               outcome  = "SCNPT",
               conditions = c("PS", "TJ", "TR", "IA", "IP"),
               incl.cut = 0.8,
               n.cut = 1,
               include = "?",
               details = TRUE, 
               show.cases = TRUE)

TS1 <- minimize(data = mydata_cpa_id,
                outcome  = "SCNPT",
                conditions = c("PS", "TJ", "TR", "IA", "IP"),
                incl.cut = 0.89,
                n.cut = 1,
                include = "?",
                details = TRUE, show.cases = TRUE)

TS2 <- minimize(data = mydata_cpa_id,
                outcome  = "SCNPT",
                conditions = c("PS", "TJ", "TR", "IA", "IP"),
                incl.cut = 0.80,
                n.cut = 2,
                include = "?",
                details = TRUE, show.cases = TRUE)

TS <- list(TS1, TS2)

RF <- rob.fit(test_sol = TS, 
              initial_sol = IS,  
              outcome = "SCNPT")
RF


# Plotting the initial solution against the test set
rob.xyplot(test_sol = TS,
           initial_sol = IS,
           outcome = "SCNPT",
           jitter = TRUE,
           all_labels = TRUE,
           fontsize = 3,
           labs = TRUE,
           area_lab = TRUE)


# Looking at case types and robustness case parameters
rob.cases(test_sol = TS,
          initial_sol = IS,
          outcome = "SCNPT")


# Obtaining cluster diagnostics
CS <- cluster(mydata_cpa_id, sol_outi, "SCNPT", unit_id = "country", 
              cluster_id = "region", wicons = FALSE)
CS


# Plotting the results of the cluster diagnostics
cluster.plot(CS,
             labs = TRUE,
             size = 8,
             angle = 6,
             wicons = FALSE)




### SET-THEORETIC MULTI-METHOD RESEARCH ###  


### Performing single within-case for sufficiency

# Typical cases
typ_out <- smmr(results = sol_outi, outcome = "SCNPT", sol = 1, match = FALSE, cases = 1, term = 1, max_pairs = 10)
typ_out

# Typical cases for each focal conjunct in a sufficient term
typ_out_fc1 <- smmr(results = sol_outi, outcome = "SCNPT", sol = 1, match = FALSE, cases = 2, term = 1, max_pairs = 10)
typ_out_fc1

typ_out_fc2 <- smmr(results = sol_outi, outcome = "SCNPT", sol = 1, match = FALSE, cases = 2, term = 2, max_pairs = 10)
typ_out_fc2

typ_out_fc3 <- smmr(results = sol_outi, outcome = "SCNPT", sol = 1, match = FALSE, cases = 2, term = 3, max_pairs = 10)
typ_out_fc3

typ_out_fc4 <- smmr(results = sol_outi, outcome = "SCNPT", sol = 1, match = FALSE, cases = 2, term = 4, max_pairs = 10)
typ_out_fc4

typ_out_fc5 <- smmr(results = sol_outi, outcome = "SCNPT", sol = 1, match = FALSE, cases = 2, term = 5, max_pairs = 10)
typ_out_fc5

# Deviant consistency cases
dcon_out <- smmr(results = sol_outi, outcome = "SCNPT", sol = 1, match = FALSE, cases = 3, max_pairs = 10)
dcon_out

# Deviant coverage cases
dcov_out <- smmr(results = sol_outi, outcome = "SCNPT", sol = 1, match = FALSE, cases = 4, max_pairs = 10)
dcov_out


### Performing comparative within-case for sufficiency

# Pairs of typical and typical cases for each focal conjuct in a term

typtyp_out1 <- smmr(results = sol_outi, outcome = "SCNPT", sol = 1, match = TRUE, cases = 1, term = 1, max_pairs = 20)
typtyp_out1

typtyp_out2 <- smmr(results = sol_outi, outcome = "SCNPT", sol = 1, match = TRUE, cases = 1, term = 2, max_pairs = 20)
typtyp_out2

typtyp_out3 <- smmr(results = sol_outi, outcome = "SCNPT", sol = 1, match = TRUE, cases = 1, term = 3, max_pairs = 20)
typtyp_out3

typtyp_out4 <- smmr(results = sol_outi, outcome = "SCNPT", sol = 1, match = TRUE, cases = 1, term = 4, max_pairs = 20)
typtyp_out4

typtyp_out5 <- smmr(results = sol_outi, outcome = "SCNPT", sol = 1, match = TRUE, cases = 1, term = 5, max_pairs = 20)
typtyp_out5


# Pairs of typical and individually irrelevant (IIR) cases for each focal conjunct in a term
typiir_out1 <- smmr(results = sol_outi, outcome = "SCNPT", sol = 1, match = TRUE, cases = 2, term = 1, max_pairs = 20)
typiir_out1

typiir_out2 <- smmr(results = sol_outi, outcome = "SCNPT", sol = 1, match = TRUE, cases = 2, term = 2, max_pairs = 20)
typiir_out2

typiir_out3 <- smmr(results = sol_outi, outcome = "SCNPT", sol = 1, match = TRUE, cases = 2, term = 3, max_pairs = 20)
typiir_out3

typiir_out4 <- smmr(results = sol_outi, outcome = "SCNPT", sol = 1, match = TRUE, cases = 2, term = 4, max_pairs = 20)
typiir_out4

typiir_out5 <- smmr(results = sol_outi, outcome = "SCNPT", sol = 1, match = TRUE, cases = 2, term = 5, max_pairs = 20)
typiir_out5


# Pairs of typical and deviant consistency cases
typdcon_out <- smmr(results = sol_outi, outcome = "SCNPT", sol = 1, match = TRUE, cases = 3, max_pairs = 10)
typdcon_out


# Pairs of deviant coverage and individually irrelevant (IIR) cases
dcoviir_out <- smmr(results = sol_outi, outcome = "SCNPT", sol = 1, match = TRUE, cases = 4, max_pairs = 10)
dcoviir_out